Objectif : Exploration des données et modélisations
Date début de l’analyse : 17 mai 2022
Date de la dernière modification : 13 juin 2022
# Packages
packages = c("tidyverse", "glue", "kableExtra", "knitr", "plotly", "hrbrthemes", "gridExtra", "stats", "arrow", "lme4", "tidymodels", "broom.mixed", "multilevelmod", "dotwhisker", "rAmCharts", "DescTools")
## Installation des packages si besoin et chargement des librairies
package.check <- lapply(packages,
FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE)
library(x, character.only = TRUE)}})
# Import BSO complet (2013-2020)
data <- read_csv("../data/4.process/managing-variables_2nd-part/whole_bso.csv", col_types = cols(author_position = "c"), na = c("", "NA"))
# Création de nouvelles variables : couleur de journal finale et APC fiables
data <- data %>% mutate(oa_color_finale = case_when(oa_color_article_BSO == "diamond" ~ "diamond",
oa_color_openalex == "gold" ~ "gold",
oa_color_openalex == "hybrid" ~ "hybridOA",
TRUE ~ "other"),
amount_apc_EUR_trusted = case_when(apc_source == "openAPC_estimation_publisher" | apc_source == "openAPC_estimation_publisher_year" ~ NA_real_,
TRUE ~ amount_apc_EUR))
Pour les modélisations, nous cherchons à expliquer et prédire le montant des APC (Y) avec 2 variables potentielles :
L’analyse exploratoire nous dira quelle variable des APC est à retenir pour inclure dans les modèles. Dans tous les cas, les variables explicatives utilisées pour prédire Y seront “oa_color_finale”, “tier”, “bso_classification”, “is_french_CA_bso_wos_openalex_single_lang”, “year” et “is_covered_by_couperin”.
Pour préparer les données nous créons donc les variables
“oa_color_finale” et “amount_apc_EUR_trusted”, puis
nous sélectionnons les 6 variables qui entrent dans les modèles en
supprimant les valeurs manquantes, et enfin nous filtrons sur les
articles pour lesquels des APC ont été payés
(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification),
apc_has_been_paid = as.factor(apc_has_been_paid),
tier = as.factor(tier),
oa_color_finale = as.factor(oa_color_finale),
# On modifie les valeurs des catégories pour éviter les conflits de nom de classes dans les modélisations
is_french_CA_bso_wos_openalex_single_lang = as.factor(case_when(is_french_CA_bso_wos_openalex_single_lang == 1 ~ "french CA",
is_french_CA_bso_wos_openalex_single_lang == 0 ~ "non french CA")),
is_covered_by_couperin = as.factor(case_when(is_covered_by_couperin == 1 ~ "covered",
is_covered_by_couperin == 0 ~ "non covered")))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == 1, amount_apc_EUR > 0)
table_y1_class <- bso_pop %>% group_by(year, bso_classification) %>%
summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup()
# Plot
graph_y1_class <- ggplot(table_y1_class, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y1_class)
table_y1_tier <- bso_pop %>% group_by(year, tier) %>%
summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(tier = as.character(tier))
# Plot
graph_y1_tier <- ggplot(table_y1_tier, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y1_tier)
table_y1_ca <- bso_pop %>%
group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>%
summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>%
mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))
# Plot
graph_y1_ca <- ggplot(table_y1_ca, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y1_ca)
table_y1_color <- bso_pop %>% group_by(year, oa_color_finale) %>%
summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))
# Plot
graph_y1_color <- ggplot(table_y1_color, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y1_color)
table_y1_coup <- bso_pop %>% group_by(year, is_covered_by_couperin) %>%
summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))
# Plot
graph_y1_coup <- ggplot(table_y1_coup, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y1_coup)
table_y2_class <- bso_pop %>% group_by(year, bso_classification) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# Plot
graph_y2_class <- ggplot(table_y2_class, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y2_class)
table_y2_tier <- bso_pop %>% group_by(year, tier) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))
# Plot
graph_y2_tier <- ggplot(table_y2_tier, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y2_tier)
table_y2_ca <- bso_pop %>%
group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>%
mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))
# Plot
graph_y2_ca <- ggplot(table_y2_ca, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y2_ca)
table_y2_color <- bso_pop %>% group_by(year, oa_color_finale) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))
# Plot
graph_y2_color <- ggplot(table_y2_color, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y2_color)
table_y2_coup <- bso_pop %>% group_by(year, is_covered_by_couperin) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))
# Plot
graph_y2_coup <- ggplot(table_y2_coup, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph_y2_coup)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable. Nous supprimons donc la variable des APC du BSO puis filtrons sur les articles ayant un APC fiable connu.
n_old <- nrow(bso_pop)
bso_pop <- bso_pop %>% select(-amount_apc_EUR) %>% filter(!is.na(amount_apc_EUR_trusted))
Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et des valeurs connues pour les autres variables.
amBoxplot(bso_pop$amount_apc_EUR_trusted, xlab=" ", ylab="Montant des APC fiables", main="Boxplot de Y")
L’étendue (écart entre les valeurs minimales et maximales) du montant des APC est importante, avec des valeurs allant de 2.27 à 9848.4 euros. Nous constatons 2 valeurs extrêmes qui se détachent largement des autres sur le boxplot, où les montants sont respectivement de 9848.4 et 9028.43 euros. Après vérification à partir de l’ISSN, les valeurs sont plausibles donc nous décidons de laisser ces valeurs aux données à modéliser.
# Normalité de Y
mean <- mean(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)
sd <- sd(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)
ggplot(bso_pop, aes(bso_pop$amount_apc_EUR_trusted)) +
geom_histogram(aes(y=..density..), color="black", fill = "steelblue", alpha = 0.2) +
geom_density( color="red", size = 1, adjust = 5) +
stat_function(fun = dnorm, colour = "Black", size = 1, args = list(mean = mean, sd = sd)) +
xlim(c(min(bso_pop$amount_apc_EUR_trusted)-1, max(bso_pop$amount_apc_EUR_trusted)+1))+ ggtitle("Distribution du montant des APC fiables") +
xlab("Montant des APC fiables") + ylab("Densité") +
theme_classic()
# Test statistique de normalité
out <- JarqueBeraTest(bso_pop$amount_apc_EUR_trusted, robust = FALSE, method = "chisq") # Y non normal (on rejette H0)
library(formattable)
df <- data.frame(value = unlist(out))
tdf <- as.data.frame(t(df))
formattable(tdf)
| statistic.X-squared | parameter.df | p.value | method | data.name | |
|---|---|---|---|---|---|
| value | 49797.6693637252 | 2 | 0 | Jarque Bera Test | bso_pop$amount_apc_EUR_trusted |

Le montant des APC fiables n’est pas normalement distribué.
# Fonction pour sauvegarder les résultats des tests
save_results <- function(output, num){
df <- data.frame(value = unlist(output))
tdf <- as.data.frame(t(df))
assign(glue("tdf_{num}"), tdf[,1:5], envir = .GlobalEnv)
}
# Test de corrélation (Spearman car distribution non normale)
output <- cor.test(bso_pop$amount_apc_EUR_trusted, bso_pop$year, method = "spearman") # corrélation significative
df <- data.frame(value = unlist(output))
tdf <- as.data.frame(t(df))
tdf <- tdf %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"), #escape special characters
data.name = gsub("bso_pop\\$", "", data.name))
formattable(tdf)
| statistic.S | p.value | estimate.rho | null.value.rho | alternative | method | data.name | |
|---|---|---|---|---|---|---|---|
| value | 242633003532012 | 0 | 0.149384487346854 | 0 | two.sided | Spearman’s rank correlation rho | amount_apc_EUR_trusted and year |
# Tests de dépendance
# Année
save_results(chisq.test(bso_pop$year, bso_pop$bso_classification), "1")
save_results(chisq.test(bso_pop$year, bso_pop$tier), "2")
save_results(chisq.test(bso_pop$year, bso_pop$oa_color_finale), "3")
save_results(chisq.test(bso_pop$year, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "4")
save_results(chisq.test(bso_pop$year, bso_pop$is_covered_by_couperin), "5")
# Discipline
save_results(chisq.test(bso_pop$bso_classification, bso_pop$tier), "6")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$oa_color_finale), "7")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "8")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$is_covered_by_couperin), "9")
# Tier
save_results(chisq.test(bso_pop$tier, bso_pop$oa_color_finale), "10")
save_results(chisq.test(bso_pop$tier, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "11")
save_results(chisq.test(bso_pop$tier, bso_pop$is_covered_by_couperin), "12")
# Couleur OA
save_results(chisq.test(bso_pop$oa_color_finale, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "13")
save_results(chisq.test(bso_pop$oa_color_finale, bso_pop$is_covered_by_couperin), "14")
# French CA
save_results(chisq.test(bso_pop$is_french_CA_bso_wos_openalex_single_lang, bso_pop$is_covered_by_couperin), "15") #toutes dépendantes
# Présentation des résultats dans une table
output <- rbind(tdf_1,tdf_2,tdf_3,tdf_4,tdf_5,tdf_6,tdf_7,tdf_8,tdf_9,tdf_10,tdf_11,tdf_12,tdf_13,tdf_14,tdf_15) %>%
mutate(p.value = as.numeric(p.value),
Dependance = case_when(p.value <= 0.05 ~ "Oui",
TRUE ~ "Non"))
# Escape special characters
output <- output %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"), #escape special characters
data.name = gsub("bso_pop\\$", "", data.name)) #remove base name
rownames(output) <- NULL #remove rownames
# Display table
kable(output, format = "html", caption = "Résultats des tests de dépendance entre les variables qualitatives") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| statistic.X-squared | parameter.df | p.value | method | data.name | Dependance |
|---|---|---|---|---|---|
| 1844.04364295786 | 63 | 0 | Pearson’s Chi-squared test | year and bso_classification | Oui |
| 10564.3561750721 | 21 | 0 | Pearson’s Chi-squared test | year and tier | Oui |
| 89.5375126042915 | 7 | 0 | Pearson’s Chi-squared test | year and oa_color_finale | Oui |
| 207.441377682062 | 7 | 0 | Pearson’s Chi-squared test | year and is_french_CA_bso_wos_openalex_single_lang | Oui |
| 155.842642000071 | 7 | 0 | Pearson’s Chi-squared test | year and is_covered_by_couperin | Oui |
| 7903.46367919093 | 27 | 0 | Pearson’s Chi-squared test | bso_classification and tier | Oui |
| 4498.17074461241 | 9 | 0 | Pearson’s Chi-squared test | bso_classification and oa_color_finale | Oui |
| 223.795592556606 | 9 | 0 | Pearson’s Chi-squared test | bso_classification and is_french_CA_bso_wos_openalex_single_lang | Oui |
| 13109.2394339774 | 9 | 0 | Pearson’s Chi-squared test | bso_classification and is_covered_by_couperin | Oui |
| 11883.1772351368 | 3 | 0 | Pearson’s Chi-squared test | tier and oa_color_finale | Oui |
| 228.437752033453 | 3 | 0 | Pearson’s Chi-squared test | tier and is_french_CA_bso_wos_openalex_single_lang | Oui |
| 15863.9058727782 | 3 | 0 | Pearson’s Chi-squared test | tier and is_covered_by_couperin | Oui |
| 1870.71182474937 | 1 | 0 | Pearson’s Chi-squared test with Yates’ continuity correction | oa_color_finale and is_french_CA_bso_wos_openalex_single_lang | Oui |
| 49839.9633692406 | 1 | 0 | Pearson’s Chi-squared test with Yates’ continuity correction | oa_color_finale and is_covered_by_couperin | Oui |
| 1145.55413773546 | 1 | 0 | Pearson’s Chi-squared test with Yates’ continuity correction | is_french_CA_bso_wos_openalex_single_lang and is_covered_by_couperin | Oui |
# Vérifier :
#- homoscédasticité (homogénéité de la variance au sein des sous-pops)
#- indépendance (sous-pops doivent avoir un effet aléatoire sur Y)
#- autocorrélation (résidus et erreurs ne doivent pas être auto-correlés)
# Transformation de la variable année
bso_pop <- bso_pop %>%
mutate(year = year - 2013)
# Modèle linéaire
lm1 <- lm(amount_apc_EUR_trusted ~ year, data = bso_pop)
# Summary du modèle
summary <- tidy(lm1)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 1633.06034 | 5.168348 | 315.97340 | 0 |
| year | 39.79084 | 1.063911 | 37.40053 | 0 |
# Y estimé
prediction_lm1 <- round(predict(lm1),2)
# Mesures d'erreurs
mae_lm1 <- mean(abs(prediction_lm1 - bso_pop$amount_apc_EUR_trusted))
mdae_lm1 <- median(abs(prediction_lm1 - bso_pop$amount_apc_EUR_trusted))
sae_lm1 <- sum(abs(prediction_lm1 - bso_pop$amount_apc_EUR_trusted))
Chaque année, le montant des APC (fiables) augmente de 39.79 euros en moyenne.
require(ggiraph)
require(ggiraphExtra)
require(plyr)
ggPredict(lm1, se=TRUE, interactive=TRUE)